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ABSTRACT 

Most workers in the field of atomic clocks encounter 
frequency and time instabilities which can be character- 
ized (or modelled) as random fluctuations. These random 
fluctuations typically display a power spectral density 
which varies as a power-law over some significant range 
of (Fourier) frequencies (e.q. , S (f) = h«f , where Y 
denotes the normalized, instantaneous frequency and f 
denotes the Fourier frequency). Typical oscillators 
and/or clocks may have regions where one specific power- 
law predominates and other regions where other power- laws 
predominate. In general, various combinations of five 
different power- laws seem to be adequate to describe 
almost all observed random behavior in atomic clocks. 

The five types are: 

White phase modulation $ (f) = h«f? 

Flicker phase modulation S'(f) = toffi 

White frequency modulation S*(f) = %f°_-. 

Flicker frequency modulation S'(f) = h_-.fi 
Random Walk frequency modulation S^(f) = h_|f L 

In addition to the random components, oscillators and 
clocks often show systematic, (i.e. , deterministic) 
trends such as offsets in frequency and time, as well as 
linear drifts in frequency. 

For the atomic clocks used in the NBS Time Scales, an 
adequate model is the superposition of white FM, random 
walk FM, and linear frequency drift for times longer than 
about one minute. The model has been tested on several 
clocks using maximum likelihood techniques for parameter 
estimation and the residuals have been "acceptably ran- 
dom." Conventional diagnostics indicate that additional 
model elements contribute no significant improvement to 
the model even at the expense of the added model 
complexity. 

I. INTRODUCTION 

Many authors (1, 2, 3) have documented the fact that most precision 
oscillators and clocks exhibit both random and systematic variations in 
their output signals. The typical random parts may include white noise 


^Deceased 


295 



phase modulation (PM), flicker noise PM, white noise frequency modulation 
(FM), flicker FM, and random walk FM. A subset of these five noises is 
usually adequate. In addition, most oscillators also exhibit a 1 i near 
drift in frequency, which is often difficult to measure. 

Experimenters often diagnose the various noises using the two-sample 
variance (or "Allan Variance") (4,5). On occasion, they will use an esti- 
mate of the power spectral density of the frequency fluctuations (4, 5). 
Of course, one cannot adequately observe the fluctuations of a single 
clock or oscillator by itself — one must look at the difference between 
two clocks. The allocation of noise levels to individual clocks requires 
three or more clocks of comparable quality. This allocation process does 
not always provide reasonable results. In fact, often the process yields 
negative values for the variance — an undesirable artifact of the estima- 
tion procedure. 

The Allan Variance is defined (1) as the infinite time average of 
sample variances based on a sample size of only two adjacent values of 

n 

1 z (y n - y n .i) 2 

R * — n - (1) 

n=l 

where y is the average frequency departure from nominal , averaged over 
the time interval and divided by the nomi nal frequency. An equivalent 
form of Eq. (1) is: 

o y 2 (x) ■ !£ ■ 1 1 W 2 (2) 

n=l 2 t 


Trequency. inac is, 


_ Lim 
°y ^ " n^> 


where x(t) and y(t) are related by 

y(t) =4 x(t) O) 

and X(t), the i nstantaneous time error, is related to the phase error of 
the oscillator by the relation: 

X(t) = (4) 

where <j»(t) is the phase error and v„ is the nominal frequency (e.g. , 
5MHz). 0 
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The Allan Vari ance is normally computed from finite data sets of the 

time difference, X„ , where 
n 

X n = X(nt 0 ) (5) 

and the estimated Allan Variance is 


1 2 / \ _ 1 

°y (T) " TF2iii 

where x = mx„. 

Although Eq. (6) is very close in form to the definition of the Allan 
Variance (see Eq. (1)), it is NOT an optimum estimator of the "true" Allan 
Variance. That is, there are other statistical techniques which provide 
more precise estimates of the frequency variability. These improved 
techniques, however, are usually valid only for very specific clock models. 
Fortunately, commercial cesium beam atomic clocks have been studied exten- 
sively, and good models are well documented. 

II. Optimum Estimates 

In the introduction, we identified two problems: 

A. Statistically inefficient estimators of the level of oscillator noises 
and drift, and 

B. Difficulties in separating individual clock performances. 

While these two problems cannot be totally eliminated, they are amenable 
to optimal estimation techniques. That is, we can minimize their effects. 

The means of estimating these parameters has been developed by R.H. 
Jones and P.V. Tryon (6, 7). Basically, the technique is that of maximum 
likelihood estimation. The technique requires an ensemble of comparable 
clocks (M > 2) and time difference data between clocks covering a signifi- 
cant duration (e.g. , a year). With the assumptions that the perturbing 
noises are both independent and Gaussian, and that the basic model is 
adequate, then it is possible to form the likelihood function as a func- 
tion of the oscillator parameters. The likelihood function is obtained 
from a Kalman Filter algorithm applied to the clock ensemble data. 

Using essentialy the same notation as used by Gelb (8), the clock model 
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and measurements can be expressed as follows: 
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where the subscripts on the matrices denote the recursion number 
(i.e. , time). 
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( 1 0 -1 0 0 0 ...) 

( ) 

(1 0 0 0-1 0 ...) 

H = ( ) 

( » . . . . ....) 

( ) 

( . . » . . ....) 


(9) 


R = 0 


(10) 


where the the number of clocks is M, the state vector, X, is a 2M column 
vector, 4> and § are 2M by 2M square matrices, and the measurement ma- 
trix, H, is M-l by 2M, since there are only M-l independent clock dif- 
ferences. 

In matrix form the equations become: 


-n 


= <t> 


-n-1 


( 11 ) 
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(12) 


and the measurements , Z„ , are: 

-n * 


2 n * H . X 

The forecasts of X and Z to step n+1 based on data up to and including 
step n are: n 

A 



= <t> 


V + 

-n 


(13) 


Ll = ‘ k+1 w 

Of i nterest are the innovations at step n+1. The innovations are given 
by ~ 

^n+1 ~ 2 n+ i " 2 n +l 
with the covariance matrix V n+ -^ 

W “ »’ • -n+1 • a 


(15) 

(16) 


where P + + is the error covariance matrix for the state vector (see 
Appendix A 1 for a brief summary of the Kalman f i 1 ter relations). 


Assuming that the driving noises, e and n n » are normal random devi- 
ates with zero mean, then the multivariate probability distribution can be 
written in the form 


f(z r Zg* • • 


The function, £, 
N 

£ = 2 In 
n=l 


) = [(27t)" l/2 |y I % I" 1 exp [-is z I ■ y' 1 - z ] (17) 

II n=0 n n n 

given by -2 times the log of the likelihood function, is 




N 

+ t 
n=l 



(18) 


2 

Now, £ is an implicit function of the parameters a, because both the 
innovations and the error covariance matrix, £~, are dependent on these 
model parameters. The estimation procedure finds that set of parameters 
(cr's) which minimizes £ (that is, maximizes the likelihood function). 
Unfortunately, £ is a non-linear function of the parameters and must be 
calculated by a complete pass through the data for each tri al set of the 
2M parameters. For example, if one has M=10 clocks and daily time differ- 
ence data for a year, then one has 365 x (M-l) = 3285 i ndependent measure- 
ments and 2M = 20 parameters to adjust in order to maximize the likelihood 
function. There exist standard computer algorithms to perform such calcu- 
lations. 
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Three additional concerns are (a) the estimates of confidence intervals 
for the parameters, (b) the diagnostics to test the adequacy of the basic 
model assumptions, and (c) the extension of the maximum likelihood esti- 
mates to include a frequency dri ft parameter for each c]ock (9). The 
model adequacy can be tested by testing of the residuals (Z ) for "white- 
ness" (i.e. , randomness); and by comparing results to morexomplex model 
assumptions. References 6,7 include a discussion of the methods used to 
estimate the confidence intervals of the parameter estimates. 

III. Experimental Results 

For many years, the National Bureau of Standards (NBS) has accumulated 
large quantities of clock comparison data on the commercial cesium clocks 
used in the NBS time scale. We used a recent sample of time comparisons 
on a dozen clocks over about two months sampled every two hours. We also 
used another set of daily data on seven clocks over a period of one year. 

The basic model assumption was that of white FM noise plus random walk 
FM noise plus linear frequency drift. Thus, for each clock in a data set 
we estimated a , a , and D the drift parameter. Also estimated were the 
corresponding Confidence intervals. The three parameters can be related 
to the more conventional Allan Variance through the equation (see Appendix 
B): 



(mc 0 ) 



a 2 (2n 2 + 1) 

+ —K~ + 

6nx 0 


(Dnx 0 ) 1 2 3 

~V 


(19) 


Figure 1 displays plots of the Allan variance obtained from the use of 
Eq. 19, above and the estimated parameters. Figure 2 displays a cumula- 
tive periodogram of residuals for one of the clocks. A periodogram of 
pure "white" noise would fall within the boundaries shown 90% of the time. 
On the shorter data run, (~ 2 mos.) linear frequency drift was not statis- 
tically significant. In fact, even on the longer run (1 year), only 
infant clocks or older clocks approaching end of life showed significant 
drift. (Of course, the algorithm could only detect relative drifts be- 
tween clocks, not a common drift shared by all clocks.) Tests were made 
using more complex models, but any improvement was found to be statistical- 
ly insignificant. 

IV. Conclusions 

A viable clock model for commercial cesium beam clocks consists of 
three elements: 


(1) White FM 

(2) Random wal k FM 

(3) Linear frequency drift 
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Maximum likelihood estimation techniques yield, reasonable results and 
confidence intervals also. Conventional tests show the model to adequate- 
ly describe observed clock behavior. Further, the technique allows one to 
estimate the individual performance of each clock. As pointed out by 
Jones, one can avoid the problem of negati ve vari ances by using a log 
transformation, y = £n(cr 2 ). 

Equation 19 al 1 ows one to express the results in the form of conven- 
tional Allan Variances. 

The new NBS time scale algorithm (TA(NBS)) makes use of the parameter 
estimation routines covered in this paper. The technique is also used for 
NBS clock calibrations. 
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APPENDIX A 


SUMMARY OF DISCRETE KALMAN EQUATIONS 


Model : 

in 


4» • JL T + U n 
-Hi -n 

Measurement: 

in 

as 

H X„ + V„ 
— — n — n 

Forecast: 

A 

X 

ss 

*k-i 

Error Covariance: 

in 

= 

$ • £j-i • + a 

Kalman Gain: 

in 

as 

£n • r • CH • ‘ 

Error Covariance: 

V 

ss 

£n ■ in " i ‘ £n 

State Update: 
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Allan Variance 
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Linear Frequency Drift 
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Allan Variance 


o (mx 0 ) = 


Mb 2 [(n+2m) 2 - 2(n+m) 2 + n 2 ] 2 

““ — iTT 2 


[!sDi 0 2 - (2m 2 )] 2 

"~^ 7 ~ 

= <s(Dmt 0 ) 2 


Composite: Assumes noises statistically independent. 

2 (2m 2 + 1) 
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°y (mt 0 ) = 


- " e + a 2 


mx 


6mx, 


+ %(Dmt 0 ) 2 


If the time error, X n , is sampled from a continuous process, then* 

2 


°y 2(m V = “*2 + 

mx^ 


a 2 mx rt 
n o 

1 7T 




* Private communication C. Greenhall. 
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QUESTIONS AND ANSWERS 


DR. STEIN: 

The model which you say seems to fit most of our clocks contains random 
walk frequency, no flicker frequency. Therefore, much much more pessi- 
mistic view of the clocks than we have usually adopted, and I was wonder- 
ing if you could comment on, for instance, what happens if you include a 
flicker term, do the residuals get better, or, what is the reason for re- 
jecting the flicker frequency model? 

MR. BARNES: 

OK. Everybody — I guess from my reputation, if anybody used flicker 
noise it would have been me, but, if you had had flicker FM present, what 
would happen in these models would be that the minimum value would move 
up and tend to have a flatter region right in the center if you really 
had flicker noise. 

I can't go and say that, unequivocably, there is no flicker noise in 
these clocks, but, I can say that if you looked at the residuals for the 
models without flicker noise and then added flicker noise, you would find 
that the improvement in the whiteness of the residuals was not statisti- 
cally significant. 


At least, that is our experience In this data. 


We tried to test the model, to see if adding other noises would 
significantly reduce the magnitude of the residual, and for no other 
model did we find significant improvement. That doesn't mean that it does 
not exist, but it does mean that in this sample, we were unable to observe 
it. 

DR. WINKLER: 

I just want to take this beautiful opportunity to point out that your fre- 
quency standard is an excellent one for long-term. 

MR. BARNES: 

I guess that's what that things says, and it doesn't have error bars on 
this particular graph, but I would guess that this is not surprising, be- 
cause the fact that it is so bad in short-term, means that it is harder 
to measure reliably the long-term performance, and hence, the confidence 
intervals along the random walk frequency, are so broad that it allows 
it to almost look too good. And that, I suspect, is an artifact of any 
kind of analysis. Thank you. 
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DR. STEIN: 

In this case, I would have to say that our qualitative observations from 
many years indicate that this analysis is probably correct. 

DR. WINKLER: 

Nevertheless, I think that it is important to note that one should not 
use short-term stability parameters for rating of clocks long-term. 

MR. ALLAN: 

Johnson tried, independently, to assess the ability for the clock en- 
semble, and the assessment has been over a year period, and a year elaps- 
ed, and another year period, so there was a great deal of difference in 
the stream of data — and the parameters on this particular clock were 
the same, two totally independent years. 
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